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The mean values of a many-body Hamiltonian including a proton-neutron pairing term and matrix 
elements of one-, two- and four-body operators within a basis of particle number projected BCS 
states, are analytically expressed in terms of a single function Q(N) depending on the number of 
particles, N. The function Q(N) is calculated using a recursion in N in which the shells and the BCS 
angles are kept the same for any step of iteration. An illustrative example is numerically considered 
in a restricted single particle space. Some specific features for the standard BCS, the projection 
after variation approach as well as for the variation after projection formalism, are pointed out. 



I. INTRODUCTION 



Short after the theory of superconductivity [l], 
showed up it has been realized [3[ that such a formal- 
ism might work also for nuclear systems although the 
constituents number is relatively small. Many applica- 
tions have been performed with a single constant for the 
interaction strength of various paired states. All calcu- 
lations were based on the supposition that the particle 
number is conserved only in the average. Amazingly the 
pairing force and the emerging seniority scheme has been 
introduced much earlier, by Racah Q. A serious ques- 
tion arose, namely to what extent the particle number 
fluctuations affect some physical observablcs. Attempt- 
ing to answer such a question many authors used a pro- 
jected particle number formalism. Two distinct calcula- 
tions have been employed. Occupation probabilities are 
determined first variationally with a standard BCS wave 
function and then the components of given number of nu- 
cleons are projected out. The resulting procedure is con- 
ventionally called the particle number projected BCS and 
as we said already the projection is performed after vari- 
ation. The second set of calculations perform the vari- 
ation after projection. The particle number projection 
has been first considered by Bayman in Ref. pj. There 
the averages of the Hamiltonian, including the mean field 
and pairing terms, and of particle number operator are 
expressed as functions of some residuum integrals which 
where estimated by the saddle-point method. The saddle 
point approximately satisfies an equation which is similar 
to the particle number equation of the BCS formalism. 
Under these conditions the Euler-Lagrange equations ob- 
tained with a particle number projected variational state 



are identical to the standard BCS equations. The pro- 
jection procedure was improved in Refs. 0, 0]. The 
residuum integrals were calculated by the saddle-point 
method with the integral path chosen such that it crosses 
the saddle-point on a line of steepest descent. Moreover 
two-dimensional recursion formulas for the residuum in- 
tegrals were provided. An extensive analysis of ordinary 
BCS, PBCS (projected BCS with projection after varia- 
tion) and FBCS (variation after projection) is performed 
within a two levels pairing model [7(. 

Another feature which was considered, referred to the 
centrifugal Coriolis anti-pairing effect. The Coriolis in- 
teraction tends to decrease the pairing strength and at a 
critical angular frequency t he g ap equation has only the 
trivial solution. In Refs. [H, [ltj it is found that crossing 
the critical point the rotational energy exhibits a discon- 
tinuity which as a matter of fact, is not confirmed by 
the FBCS calculations of Rcf. 0. Clearly the occupation 
probabilities emerging from a FBCS formalism arc differ- 
ent from those associated with the ordinary BCS, the dif- 
ference being function of the pairing interaction strength. 
The idea of particle number projection was extended to 
angular momentum. Indeed, the cranking model with 
a particle number projected wave function was consid- 
ered @ to investigate the backbending phenomena. The 
simultaneous projection of the particle number and an- 
gular momentum from a pairing correlated many body 
system has been considered in Ref. ll| for light nuclei. 
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Many papers have been devoted to the issues mentioned 
above [12| focusing on explaining some properties like, 
gap parameter, moment of inertia, spectroscopic factors, 
pairing versus nuclear deformation, and angular momen- 
tum. 

Although not as extensively, the Cooper pairs of one 
proton and one neutron have also been investigated 
18l | . The results reported there demonstrate that the 
generalized Bogoliubov-Valatin (BV) transformation in- 
cluding pp, nn and pn pairing is appropriate for treating 
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the pairing correlations in a self-consistent way, in spite of 
some earlier pessimistic views on this issue [ltl [2(| ■ Note 
that in a generalized (BV) formalism both the total num- 
ber of nucleons, the isospin third component (T3) and 
isospin (T) are not conserved. Therefore a simultaneous 
projection for all the three quantum numbers is neces- 
sary. This type of projection has been considered by sev- 
eral authors both numerically and analytically (2ll [22| . 
The intestability of a pn pairing model was treated in 
Refs. [23, [H[ by different methods. Thus the pairing 
Hamiltonian introduced by Richardson in Ref.[25( was 
considered in the context of the quantum inverse scatter- 
ing method. It is proved that the model is integrable by 
constructing explicitly the conserved commuting opera- 
tors. The eigenvalues of these operators were determined 
in terms of the Bethe ansatz and finally an expression for 
energy eigenvalues was possible [23|. A different method 
is applied, in Ref.(24|, to the same Richardson model 
which includes isospin-symmetry breaking terms. 

It seems that the common mathematical content with 
the above pairing correlated system is a serious ground 
for considering the pairing effect for other systems. Pair- 
ing forces acting among the quarks in two-color QCD 
matter lead to the color superconductivity [26|, [27}, as 
also confirmed by simulations in lattice gauge theories 
[28^ Superconductivity of metallic nanoparticles is dis- 
cussed in Refs. [29l - l31j . Due to its paramount impor- 
tance it is worth simplifying the formalism applied to 
other degrees of freedom following the successful path 
used for pairing of alike nucleons. 

In this context the present paper considers the isovec- 
tor pn pairing interaction with the projected total num- 
ber of particles. We aim at obtaining tractable equations 
for the residuum integrals and finally for the norms and 
matrix elements of projected states. 

This study is organized as follows. In Section II we 
study the factorization procedure applied to the expo- 
nent of a linear combination of the su(2) algebra gener- 
ators. The motivation for this investigation is the fact 
that the BCS function could be obtained by transform- 
ing the particle vacuum state with such operator. This 
is shown in Section III. The particle number projected 
function is described in Section IV. Analytical results of 
the average values of various interaction terms are given 
in Section V. Numerical results for a pairing Hamilto- 
nian considered in a restricted single particle space are 
given in Section VI The final conclusions are collected in 
Section VII. 



II. FACTORIZATION OF THE ROTATION 
OPERATOR 

In the theory of superconductivity bi-linear forms of 
creation and annihilation operators of fermions satisfy 
the commutation relations for the generators of rotations. 
For this reason, we shall first derive some very useful, in 
the theory of superconductivity, algebraic relations for 



the rotation operators. Although some of them are well 
known we present them for the sake of completeness. 

In quantum mechanics, the rotation of wave function 
is given by real angle 9 and the real unit vector n 



= cxp(-i6»nJ)*, 



(2.1) 



where J are the generators of rotation, given in Carte- 
sian coordinates. We are interested in various equivalent 
representations of the rotation operators. Concretely, we 
wish to establish a connection between the rotation pa- 
rameters 9, n from Eq. (|2.1|) and those denoted by 
and (3^ (/1 = ±,0) which define two independent factor- 
ized forms for the same rotation operator 

exp(-i6»nJ) = U= (2.2) 

~ exp(i(9n_ J + — i9noJo + i9n + J ) 

= exp(ia_ J + ) exp(— iaoJo) exp(ia+ J_) 

= exp(i/3 + J_) exp(— i/?o Jo) exp(z/3_ J + ). 

The indices show the components of coordinates in the 
cyclic basis: 



nJ 



-n_ J+ + no Jo — n + J_ 



(2.3) 



where 



J+ 
J_ 
Jo 



-1 
1 



(J1 + U2), 
(Ji - iJ 2 ), 



(2.4) 



Throughout the present paper the used notations are con- 
sistent with those of Ref . |32j . The spherical components 
of the unit vector are: 



n± = =F— -= sin(6) exp(±i$), 
v2 

n a = cos(6), 



(2.5) 



where and $ are the polar and azimuthal angles. The 
unit length of the vector n implies, in the cyclic basis, 



2n_n + + n = 1. 



(2.6) 

The decomposition (|2.2|) resembles the representation of 
rotation as a product of three rotations described by the 
Euler angles. The difference is that here all three rota- 
tion generators arc involved, and they are non-Hermitian 
operators. 

The relation between the parameters 9, n and a M can 
be found by taking the derivative over the 9 of the left 
and right sides in Eq. (|2.2[) and using the commutation 
relations for the rotation generators. The derivative over 
9 gives 

(m_ J + — in J Q + in + J_) exp(— i9n3) = 

ia'_J + exp(ia_ J + ) exp(— iag Jq) exp(ia + J_) 
— ia' exp(ia_ J + ) Jo exp(— iao Jq) exp(ia+ J_) 
+ia' + exp(ia_ J+) exp(— iao Jo) J- exp(ia + J_). 
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Here the symbol ""' stands for the derivative with re- 
spect to the variable 9. Using the commutation relations 
satisfied by the angular momentum components [32j 



[J+,J-] = -Jo, 

[J±,M = tJ±- 



one finds 



cxp(ia_ J+) Jo exp(— ia_J+) = Jo — ia~J+, 
cxp(— iaa Jo) J_ exp(zao Jo) = J-exp(io;o), 

exp(«'a_ J_)_) J_ exp(— ia_ J_|_) = J_ — ia_ J - 



(2.7) 



(2. 



In deriving the above relations we made use of the fol- 
lowing equation: 



e A Be- A 



[A,...[A,B] 



which holds for any operators A and B. One thus arrives 
at the equation 



- J+ — iriQ Jq + in_|_ J_ = ia'_J + — ia (Jo — ia_ J+) + ia+ exp(iao)( J_ 



-Jn- 



-J4 



(2.9) 



r 



that can be split into three equations 

o?_ 

ri- = a_ + ia a- — a' + exp^cto)-^-, 

no = a' Q + ia' + exp(iao)a_ , 

n + = a, exp(iao), (2-10) 

which can further be brought to the equivalent form 



ol_ + inoQ!_ + 

oIq + in+a.- , 
a+ exp(io;o)- 



n+az 



(2.11) 



Using the boundary conditions 

o_(0) =ao(0) =o+(0) =0, 

we obtain 



a±(fl) 



n ± 2tan(0/2) 
1 + -motan(6»/2)' 1 ' J 

-2i ln(cos(0/2) + in sin(6>/2)). (2.13) 

The result for the functions a ^(8), depends only on the 
commutation relations of J^ and therefore is the same 
for any representation of the SU{2) algebra generated 
by the angular momentum operators. Let us make use 
of this remark. Consider the fundamental representation 
and replace the operators 3 ^ by the Pauli matrices di- 
vided by 2. Next, we expand the left and right sides of 
the expression (|2.2[) in the parameters 9 and a M . The 
decomposition over a± is necessarily truncated at the 
terms linear in a± due the equation J|. = valid for the 
spin-1/2 representation. Thus, one obtains a system of 
algebraic equations for ot^O) whose solutions as can eas- 
ily be checked are given by equations (|2.12p and (|2.13l) . 



The algebraic method is fully equivalent to solving the 
system of ordinary differential equations (|2.10l) . 

We assumed, so far, that the rotation angle 9 and the 
unit vector n are real quantities. However, in the above 
development the condition for the mentioned variables 
to be real was not explicitly used. Equations (|2.12p and 
(|2.13[) . therefore, can be analytically continued to com- 
plex values of the parameter 9, the vectors n could also 
be complex under the condition that their square is equal 
to unity, i.e., n 2 = — 2n+n_ + Uq = 1. Obviously, these 
are the most general conditions. Any operator that can 
be written as the exponential of a linear combination of 
generators of the rotation with complex coefficients can 
be represented using a complex parameter 9 and a unit 
complex vector n. 

In calculating the state norms and the matrix elements 
of operators one needs to know another factorized form 
for the rotation operator where the factor operators show 
up in a reverse order as compared with that given in 
the second line of Eq. (|2.2[) . The reverse order is shown 
explicitly in the third line of Eq. (|2.2j) . We thus are 
looking for the rotation parameters /? M as functions of 9 
and n. 

In the derivation of equations (|2. 12|) and (|2 . 1 3[) we 
have used only commutation relations for the genera- 
tors of rotations. The explicit form of f}^{9) can be ob- 
tained by using the fact that under the replacement of 
(J+, J_, Jo) — s- (J-, J+, — Jo) the commutation relations 
(|2.7[) remain unchanged. From this we immediately find 



(3±(8) 



n±2tan(0/2) 
1 - t?7. o tan(0/2) : 



(2.14) 



/3 O (0) = 2iln(cos(0/2) - m o sin(0/2)). (2.15) 

If we are able to express /3 M through the a„, then we 
can change the order of the exponents with different op- 
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orators J M . These relations are, however, simple: 

P± = -a* (2-16) 



ft 



(2.17) 



Equations (|2.16j) and (|2.17j) are the necessary and suf- 
ficient in order that the operator U is unitary. They are 
therefore valid for real 9 and for real unit vectors n, that 
is for pure rotation. These equations cannot be contin- 
ued analytically to the complex values of 9 and n, because 
they involve the operation of complex conjugation. In the 
case of complex rotations, one should return to equations 
(|2.12[) and (|2.13[) . We express the complex parameters 9 
and n through a^, and substitute these expressions in 
formulas (|2.14p and (|2.15p . Simple calculations give 



ft 



a± exp(iao) 



1 + |af + a_ exp(iao) ' 



(2.18) 



ft = a + 2iln(l + -a + a_cxp(ia )). (2.19) 

These formulas establish the connection between arbi- 
trary complex parameters a M and ft. The inverse rela- 
tions have the form 



a± 



ftt exp(-ift) 



\0+0- expHft) 



(2.20) 



Q-o = ft-2iln(l + -/3+ft.exp(-ift))- (2.21) 

The factorized expression of any rotation like oper- 
ator is known in the literature as the Baker- Cambell- 
Haussdorff formula. The general necessary conditions 
which make this factorization possible are discussed in 
Ref. [13. 

Inverse relations allowing de-factorization of the rota- 
tion look as follows: 

cos(0/2) = cos(ao/2) + exp(iao/2) 
= cos(/3 /2) + cxpHft/2), 



n±sm(9/2) = -a± exp(zao/2) 

= \fi± expHft/2), 



(2.22) 



nosin(#/2) = sin(ao/2) + -a + a_ exp(iao/2) 
= sin(ft/2) + l -p+f3^ expHft/2). 



III. BCS STATE 

The BCS wave function can be written as a unitary 
transformation of the vacuum state. Such a form sig- 
nificantly simplifies the calculations and gives physical 
meaning to the algebraic transformations involved in the 
formalism. 



In the analogy with Eq. (|2.1I) one can write 

\BCS) = e~ iF \Q), (3.1) 

where 



F 



75EUE 



V2 



C am^OTh + Z a d am C am 



(3-2) 



is a Hermitian operator while |0) denotes the bare vac- 
uum state. Here, c^ m and S am are particle creation op- 
erators for protons and neutrons respectively, a is the 
index numbering shells. For example, in a spherical shell 
model a is the set of quantum numbers (njl), where n 
is the radial quantum number, j is the total angular mo- 
mentum, and I is the orbital angular momentum, rh de- 
notes the time reversal sub-state: d^ mh = ( — y~ m d l 



t 



34 1 and for a generalized pn 
221 ] analogous, otherwise dif- 



For alike nucleons in Ref. 
pairing interaction in Ref. 
ferent, unitary transformations have been used. 

The creation operators images through the transfor- 
mation exp(— iF) define the quasiparticle creation oper- 
ators: 



e~ iF J e iF = c f 
e- lF S e iF - y f 



(3.3) 



with: 



l«a|- 



cos(^=r)cT m - i— — sin(— )d 



V2 



\z a \ V2' 



0t m = cos(^)4 m - sin(^)c Q ™. (3.4) 



V2 



V2- 



We recognize here the Bogoliubov-Valatin (BV) transfor- 
mation: 



c f - U p r^ - V p d - 

'•am ~ u a^am >'a u 'am, 
— TT n rf — V n r - 



'otn — ^a4m V a c am- (3-5) 

Using the polar representation of the complex variable z, 



1 



-iip a 



(3.6) 



the occupation and non-occupation probability coeffi- 
cients, V and U respectively, are expressed as follows: 

K = cos(^), 

V r = e i ^+Vsm{^),T = p,n, (3.7) 

The complex variable z a can be interpreted as the fi = — 1 
component of a representative vector from the associated 
classical phase space. From Eq. (|3 . T[) it comes out that 
the BV transformation coefficients satisfy the normaliza- 
tion conditions: 



U% + \V a \ 2 = l, 



(3.8) 
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which assure that the quasiparticlc operators obey anti- 
commutation relations specific to fermions. Since = 
U™ and V£ = V™, in what follows we shall omit the 
isospin index for the BV transformation coefficients. 
Applying the operator exp(— iF) on the obvious equa- 
tions c Qm |0) = d am |0) = one obtains c am \BCS) = 
^>am\BCS) = 0, which expresses the fact that the BCS 
state is a vacuum state for the quasiparticlc operators. 

Using the proton and neutron creation operators one 
could define the bi-linear operators 

Ja+ = -f= ^""^ c im^im; (3-9) 
^ m 

J a — ~7= ^ C-CLindarhj (3.10) 
* m 

^oO = X y ] { C am C am ~ d am d) am } , (3-11) 



and J a o- The interpretation of the quasispin comes out 
from the following obvious equations: 



J a — I J a 7 J a ) 
Jao\Ja, Jao) 



= o, 



Jao\Ja, JaO) 

N 1 



(2j a + l) \J a ,J a0 ), 



JaO — 



N \, 

y-2^. + l). 



(3.13) 



Since the values for N range from to 2(2 j a + 1) it 
results that the minimum and maximum values of J a o are 
— (2j a + l)/2 and (2j Q + l)/2, respectively. Consequently, 
the quasispin has the expression: 



Ja 2 



!)• 



(3.14) 



which satisfy the SU (2) algebra commutation relations 
(cf. Eqs. (22)) : 



[J a +, Jfj-\ 

[J a ±, J/3o] 



— 8 a /3 J a Q, 

TS a pJ a ±- 



(3.12) 



We conventionally call these operators as quasispin op- 
erators. Indeed if we replace the d operators by c the re- 
sulting algebra defines the proton scnioritiy states. More- 
over replacing the c operators by d one obtains the neu- 
tron quasipin algebra which defines the neutron senior- 
ity scheme. Although there is a danger of mixing them 
up with angular momentum operator we use the nota- 
tion J for the quasispin algebra operators. Actually the 
pn pairing operators and the total number of nucleons 
form a representation of the SU(2) algebra which is dif- 
ferent from that generated by angular momentum com- 
ponents. Due to Eq. (|3.12[) ,thc transformation e~ %F ac- 
quires the significance of a quasirotation. Within this 
context the quasiparticlc operators appear to be the re- 
sult of a quasirotation applied to the creation and anni- 
hilation particle operators. 

It is worth mentioning some useful properties: J a ± = 
— Jo!=p, Jqo is a Hermitian operator. Due to the Pauli 
principle both J a+ and J Q _ are nilpotent: 



{Ja-) 



(Ja+) 



2j Q +2 



0. 



Acting on the vacuum state, one has 



Ja-|0) 
Jao\0) 



0, 



-0« + tt)|o>. 



For a given a the eigenstates of the proton-neutron pair- 
ing Hamiltonian in the restricted single particle space 
can be expressed in terms of the irreducible representa- 
tions of the SU(2) group. These are the states \J a , J a o) 
which are simultaneous eigenstates for the operators 



Thus, the state with a minimum quasispin projection to 
the z axis is a kind of Hartree-Fock vacuum for the low- 
ering quasispin operator. Also, the component z of qua- 
sispin is related to the state angular momentum and the 
total number of particles, TV, distributed on the substatcs 
of the single shell, j a , while the quasispin is given by the 
semi-degeneracy of the given single particle state. 

The transformation e~ tF is very useful for calculating 
operator matrix elements either in the quasiparticle or 
in the particle representation. As an example, we find 
the matrix element of a unit operator between the states 
cxp(i/3 Q _ J Q+ )|0) and (0| exp(i/3 a+ J Q _). 

By means of Eqs. (|2.20j) and (|2.21j) . we obtain 

(0|exp(i/3 Q+ J Q _)exp(i / 3 Q _J„ + )|0) = 
= (0| cxp(ia Q _ J a+ ) cxp(-ia QO <4o) exp(ia a+ J Q _)|0) 



| exp(-ia aO ^Qo)|0) 



cxp(ia a0 (j a + ^))|0) 



1 + ^Pa+Pa- 



2j a +l 



(3.15) 



In the above case, the parameter /3 Q o is equal to zero; 
we used Eq. (|2.21l) to express a a o in terms of /3 Q + and 

Under the unitary transformation, the generators J^ a 
become: 



-iF t iF 



where 



-Ja— — 



tiaO ~ 



fK / j am v am ' 

* m 

* m 



(3.16) 
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Under the action of a unitary BV transformation, the where 
operator exp(— iF) maps onto itself, since 



5 = e~ iF Fe iF 



' V2 
= F. 



a \ m 



F (f) = - E ( z ^ 2lV J a + ~ <e" 24V J*-) (4.5) 
\ and F = F(0). 

J Using the representation ()3.6|) . the operator F(tp) can 

(3.17) be written as a sum of scalar products of vectors p a n a 
and quasispin operators J a : 



The reciprocal relation for Eq. (|3.3[) takes the form 

dim = ^t m e"« (3.18) 



or, explicitly, 



($ — U -\-V c ~ 



(3.19) 



IV. PARTICLE NUMBER PROJECTION 



F{<P) = - (n a -J a+ + n Q+ J Q _) . (4.6) 



Here, n a are unit vectors defined by the cyclic coordi- 
nates 



n a ± = T—^exp(±i(ip a -2ip)), 
n a0 = 0. 



(4.7) 



The projection operator to the state with a definite 
number of particles is given by 



2tt 

Pn= f ± e i(N-N)^ 
J 2lT 



(4-1) 



One can check that the signature property of a projection 
operator is satisfied: 



P N P N = P N . 



(4.2) 



The particle number operator can be expressed in terms 
of the operators J a o 



where 



N a = 2J a0 + 2j a + 1. 



(4-3) 



The operator P/v acting on the BCS wave function 
gives a state with a definite number of particles 

\BCS,N) = C N P N \BCS) 



2ji 



J g e i(iV-iV)*, e -iF| 0) 



2tt 



Cn I ^ e -^ e -^)|o), (4.4) 

Z7T 



Using Eq. (|2.2p . we obtain 
exp(-iF((p)) = exp(i^T^p a (n a -J a+ + n a+ J a -)) 

a. 
a. 

x cxp(-z a a0 J a0 ) cxp(i ^ a a+ J a _), 

a a 

where, according to Eqs. (|2.12[) and (|2.13p . 



a a ± = n Q± 2tan(p Q /2), 
a a0 = -2zln(cos(p Q /2)). 



(4.8) 



Acting with the operator exp(—iF(tp)), after the fac- 
torizing it, on the vacuum state, we obtain 



exp(-zF(^))|0) = rj( cos2 0W 2 )) 



exp(i ^2 np-2tan{pp/2)Jp + )\<S). 



Now, we are in a position to find the projected state 
with N nucleons. 
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2tv 



\BCS,N) 



C 



f ^ exp(-iNip) J] (cos 2 (p Q /2)) jQ+1/2 cxp(z cxp(i(-Va + 2^))72 tan(p a /2) J a+ )\0) 
o 

Ci V n( cos2 (^/ 2 ))'' a+1/2 / ^ 7 ^ T cxp( l C 2 ^exp(~ l ^)V2tan(p Q /2)J a+ )|0) 



C w n( cos2 (^/ 2 )) 




y^exp(-iV> a )y / 2"tan(p a /2)J Q , + J |0). 

a / 

I 



The contour C encompasses the point £ = 0, the integra- 
tion is performed in the direction of increasing tp = arg£. 
By the condition (|4.2[) , the problem of finding the nor- 



malization constant reduces to calculating the overlap of 
unprojected and projected states: 



J 



(BCS,N\BCS,N) = C N (BCS\BCS,N) =C 2 N Y[(cos 2 (p a /2)f 3c " +1 J 



2j a +i f d( 1 



2ni ( N+1 



exp(i ^ exp(iip a )V2 tan(/j Q /2) J a -) exp(i( 2 exp(— itp a )V2 tan(p Q /2) J a+ )|0) 



c ^ / ^^n( cos2 (^/ 2 )+c 2sin2 ^/ 2 )) 



2jo+l 



(4.9) 



From the above equation we obtain 



C 



N 



d( 1 
2^C^+ 



r [](cos 2 (p«/2) + C 2 sin 2 (p Q /2)) 



(4.10) 



r 



In the derivation of the expression (|4.9p . we used the 
Eqs. (|2.20[) and f|2 . 2 1[) . which allow to change the order of 
the factors of the unitary operator. Also, we had to calcu- 
late the average of the product of exponents of operators 
associated to individual shells. If the two involved shells 
are distinct the exponents commute with each other and 
their average over the vacuum is unity. An average value 
different from unity occurs only when a' = a" , i.e. only 
when the two shells coincide. The average value splits 
thereby into the product of the average values for the 
individual shells. 

Further simplification can be achieved by using a bino- 
mial formula and then evaluating the integrand residue 
in the point £ = 0. After applying the binomial for- 
mula for the sum of a large number of terms and then 
finding the residue of the integrand the problem becomes 
combinatorial in nature, which is not attractive from the 
computational point of view, because the number of op- 
tions needed to be considered and the number of terms 
in the sum grow with ./V exponentially. 



From the computational point of view, the possibility 
of reducing the problem to evaluation of a recursion looks 
more attractive. We introduce the function 



Q{N) = C 



N 



with C N defined above. 



Integrating by parts, one finds 



Q(N) 



^ J 

P c 



(4.11) 



(4.12) 



where 



Qpi.0 = 



C 2 sin 2 0V2) 



N cos 2 ( P/3 /2) + C 2 sin 2 (^/2) 
[](cos 2 ( jOQ /2) + C 2 sm 2 ( Pa /2)) 



2jo + l 



(4.13) 



with 0« = 2(2 jp + 1). We expand further the expression 
in front of the product sign in powers of £ 2 . Each member 
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of the series is a function of Q(N') for some value of 
N' < N. The function Q(N), therefore, is expressed 
as the sum of Q(N') evaluated for a smaller number of 
particles. It only remains to fix the boundary value for 
N = 0. From the definition of Q{N) one easily finds 

0(0) = n(cos 2 (p Q /2)) 2 ^ +1 . (4.14) 

a 

Note that for negative integer values of N, the function 
Q(N) is equal to zero, as can be seen from the Cauchy 
theorem related to the contour integral. Also, Q(N) = 
for N = 1 mod(2). We thus get a recursion 

Q(N) = ^Q P (N), (4.15) 

o N/2 

Q P (N) = ^^(-)" +1 tan 2 "(^/2)g(7V-2r l ). 

n=l 

(4.16) 

The number of operations to calculate Q(N) grows with 
increasing N only quadratically. From the viewpoint of 
numerical calculation, estimates for oscillatory contour 
integrals are associated with considerable difficulties. We 
avoid this difficulty by reducing the problem to the com- 
putation of the recursion relations. Within a variation af- 
ter projection procedure the angles pp are determined by 
the equations provided by the conditions that the energy 
of the system for a fixed N be minimum. The same angles 
are however used to calculate the factors Q(N — 2n — 2) 
involved in the summation operation over n. The BCS 
angles as well as the number of involved shells are pre- 
served during the iteration process. 

The product factors from Eq. (|4.13p can be expanded 
in power series of l/£. In this case we obtain a recur- 
sion for calculating the function Q(N) starting from large 
numbers of particles. Recursion of this form is more con- 
venient to calculate Q(N) for values of N close to the 
maximum 

= ^^/3. (4.17) 

& 

Using the the expansion in l/£, we obtain Eq. (|4.15[) with 
Q>(N) = g2L (4.18) 

(fi-iV)/2 

x Yl (-r +1 cot 2n ( P/3 /2)Q(iV + 2n). 

n=l 

These expressions should be supplemented by the bound- 
ary condition 

- n( sin2 (-°«/ 2 )) 2jQ+i - ( 4 - 19 ) 

a 

For N > CI, the function Q(N) is identically zero. 



The only singularity of the integrand in the expression 
(|4.12j) is a pole at C = 0. For this reason, we can deform 
the contour of integration, squeezing it around zero or 
moving it to infinity. In the first case, the expansion in 
powers of £ is the appropriate one, while in the second 
case, the expansion in powers of 1/C is valid. Obviously, 
the results (|4.16[) and (|4.19|) coincide. Also we note that 
for N = the function QP(N) given by Eq. (|4.16[) is not 
defined. For this case one should use the 1/C expansion 
from Eq. (|4.19l) . By contrast for calculating Q^{Cl) the 
expression 14.161 is the appropriate one. From the inte- 
gral representation of Q(N) as well as from Eq. (A. 11) 
considered for the diagonal case, a very simple relation 
follows: 

n 

E Q( N ) = L (4- 2 °) 

7V=0 

In virtue of this expression Q(N) acquires the significance 
of the admixture probability of the iV-projected state in 
the BCS wave function. 

The average values of any fermion operator corre- 
sponding to the particle number projected BCS state, 
can be expressed in terms of the function Q(N), which 
is the inverse of the normalization constant squared. In 
next Section this will be shown for few examples. 



V. AVERAGE ENERGY 

The projected state might be used as a variational state 
for a proton-neutron pairing Hamiltonian. We shall cal- 
culate the average of a many body Hamiltonian with two 
body interaction whit a strength depending on shells 

# = \)N a + Y^V a pJ a+ Jp-. (5.1) 

Within the standard BCS formalism, A is the Fermi sea 
level (chemical potential) which is to be fixed by solving 
the BCS equations. Here it is just a parameter involved 
in the gap energy equation. The two terms entering the 
microscopic Hamiltonian will be separately treated. The 
situation when V a p = —2G can be used to describe a pos- 
sible transition of the pn system from the normal to the 
superconducting phase. Also such a Hamiltonian could 
be used for describing the rate of double beta Fermi-type 
decay. Note that such a particular form of the pn pair- 
ing interaction is not invariant under the rotations in the 
isospin space, but preserves the third component, T3, of 
the total isospin. Indeed, the two body interaction com- 
prises terms of isospin 0,1 and 2. Therefore the eigen- 
states of H are expected to be a mixture of components 
of different isospin. In order to have an isospin invariant 
Hamiltonian we have to account also for the pp and nn 
interaction. Due to these features we consider (|5.ip as an 
illustrative example which allows us to describe the main 
ingredients of the present formalism. 
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A. Mean field term 

The mean field energy is determined by averaging the 
particle number operator for each shell. 



(N a ) = (BCS,N\N a \BCS,N). 



(5.2) 



With the interchange of the order of the exponent oper- 
ators, as described in the previous Section, this average 
is transformed into 



where 

7VC) = 2*n(cos 2 (p 7 /2)) 2 ^ +1 

7 

. ^~(0| exp(z ^ exp(i?/' 7 )v / 2 tan(/? 7 /2) J 7 _) 



(5.4) 
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x exp (-ix (J aQ + j a + 1/2)) 

x exp(iC 2 ^cxp(-# 7 )\/2tan(p 7 /2)J 7+ )|0) L=o- 

7 

The derivative over x is taken at x — 0. Making use of 
Eqs. (|2~3|) . (j2~2"0j) and (|2~2T|) . V a {() can be simplified to 
give: 



(5.5) 



The function Ga{() enters the definition of Q(N) and is 
given by Eq. (|4.13[) . Combining Eqs. (|5 .3[) and (15. 5[) . one 
obtains 



(iV a ) = AQ^AOQ" 1 ^). 



(5.6) 



The sum over a in Eq. (|5.6p gives identity N — N. Using 
Eq. (|4.3p . one can find the average of J a p. 



B. Energy gap function 



previous subsection, we present results for the matrix 
elements diagonal in the BCS angles: 



1 



(J a+ ) = -iQ- L/A (N + 2)Q-^'(N)^=(N + 2) 

x cxp(# a )cot( j o Q /2)g Q (Ar + 2). (5.8) 

Alternatively one may express the above matrix elements 
as a polynomial in tan(p/2) by using a power expansion 
in C for the integrand, as explained already before. The 
result is: 

(J a+ ) = iQ~ 1 ' 2 (N + 2)Q~ 1 ' 2 (N)(2 ]a + l)^= (5.9) 

N 

x expHVa) t&n 2n+1 (p a /2)Q{N - 2n). 

71=0 

The matrix element of J Q _ can easily be found by com- 
plex conjugation 

(J a _) = (BCS,N\J a -\BCS,N + 2) 

= -(J a +)*- (5.10) 

This quantity multiplied by y/2 defines the spectroscopic 
factor for a pair a of states which could be measured 
in a deuteron transfer reaction. If the two body interac- 
tion strength is state independent and equal to — G/4 the 
pairing interaction term resembles the pairing interaction 
for alike nucleons. In virtue of the particle number con- 
servation, the average of GV%^2 a Ja+ is equal to zero. 
By analogy with the case of alike nucleon pairing, we call 
the sum 

^ = GY,(BCS,N + 2\cl n J afh \BCS > N), (5.11) 



the gap parameter for the A^-nuclcon system, which 
might be a good definition at least in the limit of large 
N. 



Within the BCS theory the expression of the gap 
energy function is obtained by averaging the operator 
J a + on the unprojected BCS state. Instead, for 
the particle number projected BCS formalism, the ma- 
trix element of the mentioned operator between the states 
with N and N + 2 particles is to be calculated. We start 
by calculating such a matrix element for each term under 
summation 

(J a+ ) = (BCS,N + 2\J a+ \BCS,N) 

= C N+2 (BCS\J a+ \BCS,N). (5.7) 

The BCS wave functions of the initial and final states can 
be different (p a ^ p' a , ip a ^ ip' a ). This case is considered 
in Appendix A. Here, following the path described in the 



C. Interaction energy 

Now, consider the proton-neutron (pn) pairing inter- 
action with the generic term: 



(J a +J-p) = (BCS,N\J a+ Jp-\BCS,N). (5.12) 

First, we lift the operators J a + and Jp- to the arguments 
of exponents by introducing two derivatives over x and 
y, evaluated in origin: 
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where 

^(C)=-n( C0S W2)) 2 ^ +1 (5.13) 

7 

x exp(i exp(ii/' 7 )v / 2 tan(/5 7 /2) J 7 _) 

x exp(— ixJ a +) exp(— J-/3) 

x exp« 2 ^exp(-^' (5 )V2tan( / 9a/2)J ( 5 + )|0) | x=y=0 . 

I 



(Ja+J/3-) 



In the standard BCS formalism, according to Eq. (3.6) 
the occupation probabilities for protons and for neutrons 
associated to the single particle state \a) are equal to 
each other. 

\V£\ 2 = ICI 2 - (5-16) 

Within a particle number projected formalism, the occu- 
pation probability is different from that defined above. 
Indeed, here the occupation probability for the pair of a 
states is given by: 

-2(J Q+ J Q _) ^sin 2 (p Q /2). (5.17) 



VI. ONE AND TWO SHELLS CASES 

In this Section we shall focus on the pairing Hamilto- 
nian 

H = y ^ {e a — A) {c^ am c am + d) am d am ) 

am 

- G C am d am d PnCp n . (6.1) 

a,/3;m,n 

In the special situation when the single particle space oc- 
cupied by protons and neutrons is restricted to a single 
j, the energy for N nucleons calculated as the average 
of H corresponding to A = with the iV-projccted state 
defined before, is a quantity growing quadratically with 
N and depending neither on p nor on ip. In this respect, 
one may say that the TV-projected state does not exhibit 
a superconducting character. However, the spectroscopic 
factors defined by Eq. (5.10) can be calculated if p and 
■0 are determined within the standard BCS. In that case, 
the spectroscopic factor is readily obtained once the func- 
tions Q(N) are calculated. 



The product of four exponents is then reordered, using 
the results of the Sect. II. 



The final expressions for a 7^ j3 and a — /3 become 



After some algebraic manipulations one finds a 
Bernoulli distribution 

= ( 2j N + /2 ) P N/2 (1-P) 2j+1 - N/2 (6.2) 

with probability p = sin 2 (/?/2), that can be written, in 
fact, immediately starting from the integral form (|4.10[) 
or using combinatorial arguments to calculate the norm 
of the particle number projected BCS states. 

Let us consider now that a number of nucleons, N, 
are distributed among two single particle states whose 
quantum numbers are specified through their angular 
momenta j\ and j'2 respectively. The potential entering 
Eq. (|5.ip . now has the form 

V aP = 2G. (6.3) 

The corresponding energies of the shells are denoted by 
ei and £2, respectively. 

A. Standard BCS 

If one neglect the renormalization of the single particle 
energies due to the pairing interaction, the energy of the 
constrained system of TV nucleons has the expression: 

E'(N) = £ 2(2j Q + l)(e« - A)K| 2 - (6.4) 

a 

with 

A = !^2(2j a + l)E/ a ^ a) (6.5) 

a 

N = ^2(2j Q + l)|y Q | 2 . (6.6) 

Of 



= Q-Wjxm + WW - (2j a + ij^TO ^^I^r^?^ 72 ^ (5 ' 14) 

4 tan (pp/2) — tan (p a /2) 

= -Q-\N)(-Q a (N) + (j a + 1/2) tan 2 "( /9Q /2)(2nj Q - l)Q(N - 2n)). (5.15) 

n=l 

I 
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The condition of minimum energy as a function of the 
BCS angles can be used to express the BCS angles in 
terms of the energy gap A and the parameter A: 



V a \ 2 

u 2 



1(1 e a - A 

2 l T ^~w+W 



)• (6-7) 



The energy gap and A can be found from the self- 
consistency condition (|6.5|) and the particle number con- 
straint (16.61): 



G 



(2 Ja + l)(e a -A) 



2j a + 1 



E 



= £(2j Q + 1) - N, 



Y- 



= 1. 



The variational problem also constrains the phases. 
Denoting by if the phase of the gap parameter 



A = \A\e iv , 



(6.9) 



one successively find arg V a = arg A or, using the relation 



7T 

2 



(6.10) 



The average energy (|6.4[) is independent of the phase fac- 
tor of V a . The absolute scale of the phases, therefore, is 
not determined. 

If the BCS equations admit nontrivial solutions, the 
system is, by definition, in a superconducting phase, its 
energy being calculated by Eq. (|6.4j) . 

In the case of a single shell, solution of the above equa- 
tions has the form 



ei-A 
A 2 

r 2 

,/ 

E'(N) 



l -G(2j + 1-JV), 



-G 2 N(4j 
N 

2(27+1)' 
~-N 2 . 



N), 



The quasiparticle energy is equal to y(2j + 1) while the 
system energy is obtained by subtracting from E'(N) the 
contribution of the constraint term: 



1 



E gs (N) = ei N - -GN(4j + 2-N) 



.11) 



In the BCS theory, therefore, the superconducting state 
exists for any number of particles. This conclusion, how- 
ever, is not supported in the particle number projected 
BCS theory. It is not difficult to see that the average 
energy corresponding to the ./V-projected BCS state of is 
independent of the BCS angles and phases, and equal to 



1 



E(N) = exN - -GN{Aj + 4 - N). 



(6.12) 



This result can be obtained with the general formulas of 
the previous Sections, where A is set equal to zero but also 
by exploiting the fact that after projecting the nucleon 
total number, in the BCS wave function survives only 
the component - (J+) w / 2 |0). The average value of the 
Hamiltonian H, corresponding to A = 0, for this compo- 
nent provides (|6.12[) . Note that the the unprojected BCS 
state is higher in energy than and the A^-projectcd BCS 
state. 

Moreover system defined by the Hamiltonian (|6.1|) re- 
stricted to a single shell is exactly solvable, since the 
components of the Hamiltonian are expressed through 
the Casimir operator of the quasi-rotation group, J 2 , and 
the quasi-spin projection on z-axis, J a o ■ In our case, the 
expression of (|6.12[) appears to be just the eigenvalue of 
the pairing Hamiltonian, while the TV-projected state is 
the corresponding cigenfunction. 



B. TV-projected BCS 

Here we consider the case of two shells calculations. 
For what follows it is useful to write the expression of 
the Q function in a more compact form. The reciprocal 
norm squared can be written as: 



0(2) 
0(4) 

0(6) 

0(8) 

0(io) 

where 

0(o) 



T 2 Q(0), 

i(T 2 -T 4 )Q(0), 
1 



(6.13) 



24 



;T| - 3T 2 T 4 + 2T 6 ) Q(0), 
(T 2 4 - 6T|T 4 + 8T 2 T 6 + 3T 2 - 6T 8 ) Q(0), 



1 5 

-^T 2fe Q(10-2fc) 



(6.14) 



k=l 



(cos 2 ( Pl /2)) 2 ^+ 1 (cos 2 (p 2 /2)) 2 ^+ 1 , 

(2. n + 1) tan"(pi/2) + (2j 2 + 1) tan' l (p 2 /2). 



The terms Q a (N) specified by (|4.16[) acquire also com- 



(6.15) 



pact forms: 






0"(2) = 
Q (4) = 


F a2 Q{0), 

^Q(0)(T 2 F a2 -F ai ), 




Q Q (6) = 


^O(0)(i(T 2 -T 4 )^ 2 - 


- T 2 F a 4 - 


Q (8) = 


Jo(0)(i(T 2 3 -3T 2 T 4 + 


2Te)F a 2 




2 (-^2 ~ Ti)F a i + T 2 F a Q 


— F aS ), 


Q (io) = 


5 

E^(-i) p+1 o(io- 
P =i 


2p). 



(6.16) 
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FIG. 1: The energy of a system with N nucleons (N = 2, 4, 6, 8 
), obtained within the BCS formalism, is plotted as function 
of the pairing interaction strength, G. 





G=0.1 MeV 


G=1.0 MeV 


N=2 N=8 


N=2 N=8 


Ebcs 
A 

EpbCS 

Efbcs 

A (N) 


1.37573 6.93772 
0.30183 0.51425 
1.25987 6.48216 
1.25969 6.48213 
0.4292 0.5753 


-8.35069 -21.38316 
3.31395 5.65241 
-9.35235 -25.38779 
-9.35235 -25.38779 
4.6 866 6.3196 



TABLE I: The ground state energy provided by the BCS, 
PBCS and FBCS formalism, respectively , given in MeV, are 
listed for two values of the nucleon total number and two 
values of the pairing strength, G. 



where 



(2 Jct + l)tan"(p Q /2). 



Using these partial results in connection with the matrix 
elements (5.6), (5.13) and (5.14) the system energy calcu- 
lated as the average value of the pairing Hamiltonian H 
corresponding to the TV-projected BCS state, is readily 
obtained. 

Making use of the results obtained so far one can cal- 
culate the ground state energies as well as the energy 
gap. Calculations were successively performed for the 
standard BCS, the projection after variation approach, 
PBCS, and the projection before variation formalism , 
FBCS. Also, the exact eigenvalues of H have been ob- 
tained by diagonalization. The input data in our calcu- 
lation are: 



FIG. 2: The energy of a system with N nucleons (iV = 2,4,6,8 
), obtained within the PBCS formalism, is plotted as function 
of the pairing interaction strength, G. Energies predicted by 
the PBCS, the FBCS and diagonalization (Exact) are iden- 
tical. This is indicated by assigning to the given curves ah 
three labels: PBCS,FBCS, Exact. 
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FIG. 3: The gap energy for a system with N nucleons 
(iV = 2, 4, 6, 8) is plotted as function of the pairing interaction 
strength, G, within the standard BCS formalism. 
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FIG. 4: The gap energy for a system with N micleons 
(N = 2, 4, 6, 8) is plotted as function of the pairing interaction 
strength, G, within the standard FBCS formalism. The gap 
function is defined by Eq. (|5.11[) . 



3 . _ 7 

3i 2' 32 ~ 2' 

ei = 1 McV, e 2 = 1.5 McV. (6.17) 

For a given N (= 2, 4, 6, 8) we solved the BCS equations 
(|6.8[) and then calculated the ground state energy (|6.4[) . 
Results for energies are plotted, in Fig. 1, as function of 
the pairing interaction strength. The uncorrelated sys- 
tem has the energy 

Enormal = Net. (6.18) 

Alternatively, with the parameters p\ , pi , ipi , ^2 deter- 
mined by the standard BCS approach, the system en- 
ergy was calculated with the general expression provided 
by the present formalism with the A-projected state. 
We conventionally call this approach as projected BCS 
(PBCS). The resulting energies are represented in Fig. 2 
as functions of G. 

For each N, further, the energy E(N) given by 
the TV-projection formalism was minimized with respect 
to the parameters involved in the trial function, i.e. 
Pi, P2, "01 j V>2- It turns out that the ground state ener- 
gies provided by PBCS, FBCS and diagonalization are 
equal to each others for any G in the range of 0.1-1.0 
MeV. The common values were used in Fig. 2 showing 
the energy dependence on the pairing strength parame- 
ter. 



The two figures mentioned before exhibit some com- 
mon features. Energies are decreasing functions of G. 
They show a linear dependence on G with the slope de- 
pending on the particle total number, N . For low values 
of G, energies are increasing functions of N, while for 
large pairing strength they are decreasing with A. In 
the four approaches, the transition from one energy or- 
dering to another is taking place for different values of G. 
The largest critical G is met for the standard BCS. To 
get a better view on the quantitative energy split due to 
its A-dependence we collected the energy values for the 
extreme values of AT, in Table I. From there one can see 
that the energy ordering for a given G is not depending 
on A. Indeed, for G = 0.1 MeV the ordering for A = 2 
and A = 8 are: 

Enormal > EbCS > EpBCS > EpBCS, N = 2,8, 

while for G — 1.0 MeV the energy order is as: 

Enormal > EbCS > EpBCS = EpBCS, N = 2,8. 

It is remarkable the fact that energies provided by the 
PBCS and FBCS approaches are practically the same. 
Only very tiny differences are noticed for G=0.1 McV, 
which might be caused by the numerical procedures 
adopted for the two cases. 

Due to the above mentioned ordering for the ground 
state energies provided by different formalism we may 
conclude that for all cases the superconducting phase 
is achieved for any G > 0.1. Another observable con- 
sidered in our study is the energy gap obtained within 
the standard BCS and the FBCS formalism, respectively. 
The results obtained for a fixed N were plotted as func- 
tion of G in figure 3 and 4, respectively. These figures 
show a linear dependence of both gaps on G. The split 
of gaps due to their N dependence is larger for A than 
for A^). Indeed, according to the data listed in Table I, 
A (AT = 8) - A (AT = 2) is equal to 0.212 MeV and 2.338 
McV for G equal to 0.1 MeV and 1.0 McV respectively, 
while A® - A<2) for G = 0.1; l.McV, amounts to 0.146 
MeV and 1.633 MeV, respectively. Note that for a given 
set of N and G we have A > A^ N \ Since the energy gap 
might be looked at as a measure of superconductivity one 
may expect that the superconductivity effects are more 
pronounced in the PBCS and FBCS than in the BCS. 

The Hamiltonian H can be diagonalized in a basis of 
definite number of particles, N 

ia 1; a 2 ) = c NlN M4) Ni/2 (44) N2/2 \o), 

ATi + N 2 = N. (6.19) 

with Cm x n 2 standing for the normalization constant. The 
matrix elements of H in this basis have simple expres- 
sions: 
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(JVi, N 2 \H\N U N 2 ) = exNx + e 2 N 2 - - (N^h + 4 - JV X ) + N 2 {4j 2 + 4 - N 2 j) , 
(iVi + 2,iV 2 - 2|fT|JVi, JV a ) = -- [(JVi + 2)iV 2 (4ji + 2 - iVi)(4j a + 4 - 7V 2 )] 1/2 , 
(JVi -2,N 2 + 2\H\N 1 ,N 2 ) = -j [JVi(JV 2 + 2)(4 7 ' 1 + 4 - ATi)(4j 2 + 2 - 7V 2 )] 1/2 . (6.20) 

I 



For each N and a fixed G, we diagonalized the above 
matrix and depicted the lowest eigenvalues. Further, 
these were compared with the energies provided by the 
FBCS and the PBCS formalisms. In this way we found 
out that the three sets of energies are identical. This is a 
nice example when the solution of two variational prin- 
ciple equations reproduce the exact ground state energy 
However, in a realistic single particle space and, more- 
over, when an isospin invariant Hamiltonian is instead 
considered, this feature does not necessarily show up. 

VII. CONCLUSIONS 

Many interesting properties of the pn pairing have been 
derived using for the BCS function the form (3.1). After 
making use of the factorization described in Section I, 
one finds out that the function (3.1) for a fixed a rep- 
resents, the coherent state of the SU(2) group generated 
by Jafi- Our attention was focused on the nucleon num- 
ber projected BCS. Since any matrix element can be ex- 
pressed in terms of the norms of the involved states, we 
started with the norm calculation. One of the main re- 
sults of the present paper is the recursive formula for 
these norms. By this equation the norm of a projected 
state with N particles is related with the norms of the 
iV'-projected states with N' < N. For pairing interac- 
tion of alike nucleons, a similar recursion formula was 
obtained in Ref. (0). The difference between the two 
recursion formulas consists of that there the recursion is 
operating in two dimensions, there are two indices which 
are iterated, while here only one index is involved in the 
recursion. 

To prove the usefulness of the obtained recursion for- 
mula several matrix elements have been evaluated. The 
one for the proton-neutron pairing operator of a given 
shell, is interpreted as spectroscopic factor for a deuteron 
transfer reaction. Being guided by the analogy with the 
alike nucleon pairing we defined a quantity which might 
be a measure for the energy gap in the particle number 
projected picture. Also, we calculated the occupation 
probability for a given state with a proton-neutron pair. 
In Appendix A the matrix element involved in the width 
of an a decay process is analytically expressed. Also, the 
matrix element for the two body proton-neutron inter- 
action in the particle-particle channel between the states 
associated with the mother and daughter nuclei involved 
in a double beta Fermi decay, are obtained in a compact 



form. 

Using a Hamiltonian including a mean field term and a 
proton-neutron pairing interaction we calculated the sys- 
tem energy as function of the particle number N and the 
parameters p a ,ipa defining the unprojected BCS wave 
function. Since many features of the paired system can 
be found also in a restricted single particle space we dis- 
cussed the simple cases of one and two single particle 
states. Since the energy associated to a single j and 
the iV-projected picture, is constant with respect to the 
BCS parameters we concluded that for this case a super- 
conducting phase cannot be reached. However, the case 
of two j is suitable for studying the pairing properties 
both in the standard BCS and projected BCS formalism. 
In the later situation we considered both cases when the 
variation is performed before and after projection. In the 
two level situation we have proved that the ground state 
energy provided by the PBCS, FBCS are the same and 
moreover equal to the exact ground state energy obtained 
by diagonalization. It is an open question whether this 
feature is caused by the restricted single particle space or 
by the particular choice o the model many-body Hamil- 
tonian. 

One of the evident limitations of our formalism con- 
sists in the fact that the trial BCS unprojected function 
allows us to describe only nuclei with N=Z. This feature 
can be however improved by adding to (3.1) two factors 
accounting for the proton-proton and neutron-neutron 
pairing, respectively. 

As an imminent project for the near future we also 
mention the extension of our formalism to isospin pre- 
serving Hamiltonians. In the second step of the formal- 
ism development we shall attempt to include in our study 
the proton-neutron T — pairing as well as the isospin 
projection. This plan is, in fact, a reflection of our be- 
lief that a new and powerful technical result might be 
decisive in unveiling new properties of the paired proton- 
neutron system. 
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Appendix A: Matrix elements of different final and 
initial particle number projected BCS states 

In various applications matrix elements of operators 
between two different BCS states are frequently encoun- 
tered. In the case of projected BCS, such matrix ele- 
ments are calculated in the manner described in Section 
V. Here we derive formulas for (BCSf,Nf\0\BCSi,Ni) 
for (9 = 1, J± a , J± a J±f3 and J^ a J±p with different BCS 
angles of the final and initial states. 

We start from the diagonal in particle number overlap 

(BCS f ,N\BOSi,N) = QJ 1/2 (N)(BCS f \BCS l ,N) 
= QJ 1/2 (N)Q; 1/2 (N)Q fl (N), 



where Qf(N) and Qi(N) are the diagonal Q-functions 
defined for the final- and initial-state BCS angles (p' a ,ip' a ) 
and {p a ,ip a ), respectively. Calculations similar to those 
carried out in Section V give give the off-diagonal Q- 
function 



Qfl{N) = / 2Vi II ( cos (^/2) cos(p Q /2) + C 2 exp(#^ - sin(p' Q /2) sin(p Q /2)) 2j ' Q+1 . (A.l) 
c a 

Integrating by parts and expanding into a series in £ or l/C, we obtain the recursion for the calculation of this function: 

Oi 

n N/2 

Q%{N) = ^^(-)"+ 1 (exp^-^a)tan^/2)tan(p Q /2))"Q /l (7V-2n) (A.2) 

n— 1 

(n-AQ/2 

= E (-)" +1 (expHC+ l Va)cot(p' Q /2)cot(p Q /2))"g /l (iV + 2n). (A3) 

I 



The boundary conditions follow from (|A. 1|) 

QfM = n( cos (^/ 2 ) cos (^/ 2 )) 2jQ+i > 

Qfiify = n( ex P(^-#a)sin(^/2)sin(p Q /2)) 2j ° +1 

a 

As in the case of the function Q(N), Qfi(N) vanishes 
outside the interval (0, f2) and for odd N. In the limit 
where the angles in the initial and final BCS states co- 
incide, we recover the result for Q(N). In particular, 
Q ff (N) = Q f (N) and Q tl (N) = Q t (N). Q a fl (N) van- 
ishes outside the interval (2, Q — 2). 



Now, consider the 2-fermion gap function 

(J a+ ) fl = (BCS f ,N + 2\J a+ \BCS i ,N) 

= Q] 1/2 (N + 2)(BCS f \J a+ \BCSi,N). 

The calculation is performed using the representation of 
J Q + in the exponential form 

J a+ = —i-jj^ exp(*Ja+a)la; = o 

and Eq. (J3TT5J) . A similar representation is used for the 
other quasi-spin components. After some algebraic ma- 
nipulations, we obtain 



(J a +)fi = ~iQ] 1/2 {N + 2)Q7 1/2 (N) — (A + 2)exp(#a)cot(p a /2)Q^(A + 2), (A.4) 

(J a0 )fi = Q- f 1/2 {N)Q- 1/2 {N){^Q a fl {N) ~ ( ]a + l/2)Q ft (N))- (A.5) 

I 

Summing the last equation over a one arrives at: which is consistent with the definition for particles num- 

2^(Joa)/i = Q] 1/2 (N)Q- 1/2 (N)Q fl (N)(N - ifi), 
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ber operator, Eq. (|4.3p . The expression for (J a -)ji is 
obtained by complex conjugation (|A.4|) . 
Consider now the 4-fcrmion gap function 

(J a+ Jp+) fl = (BCS f ,N + 4\J a+ Jp + \BCSi,N). 

This matrix element represents the width of the a de- 
cay of the mother nucleus in the state \BCSf,N + 1) 
to the daughter nucleus in the state \BCSi, Nf) . Also, 
based on this matrix elements one may define the spectro- 



scopic factor of a reaction which removes an a particle 
from the states a and /3. In other words squaring this 
matrix element one obtains the occupation probabilities 
of the states a and ft with proton-neutron pairs. The 
calculation uses the representation 

J a +Ji3+ = ~ exp(i,J a+ x + iJp+y)\ x=y=0 

and the method of Section II. For a ^ (3, we obtain 



(J*+j0+)fi = -Q f 1/2 (N + 4)Q- 1/2 (N)±(N + 2)((2j> + 1)Q%(N + 2) - (2j a + 1)Q%(N + 2)) 

exp(#^ + itf a ) tan(^/2) tan(^/2) 

exp(iip' a - iip a ) tan(p^/2) tan( / o Q /2) - exp(iip' - iipp) tan(p^/2) tan(p / j/2) ' 

Similar expression is obtained for the matrix element 



(A.6) 



(J a+ Jp-) fi = {BCS f ,N\J a+ Jp-\BCSi,N) 

1 

exp(iip'a — iip a ) tan(/°fl/2) tan(p a /2) 



Q] 1/2 (N)Qi 1/2 (N)jN((2j + \)Q%N) - (2j a + \)Q%{N)) 



exp(iip' a - iip a ) tan(p^/2) tan(p Q /2) - exp(i^ - iipp) tan(^/2) tan(p /3 /2) ' 



(A.7) 



We may ask oursclf how this situation may appear? The answer is offered by the double beta decay with/without 
neutrinos in the final state. Indeed, in such a process the nucleus \N,Z) goes to the nucleus \N — 2,Z + 2), two 
electrons and two/zero anti-neutrinos. The mentioned states are described by different sets of BCS angles and phases 
but have equal number of nucleons. As a matter of fact the two body interaction whose matrix element is calculated 
is nothing else but the particle-particle interaction of the Fermi type. These comments prove that to calculate such 
a matrix element is an important step in describing some important physical processes. 

The average (|A.6j) is symmetric under interchange of the indices and (J Q +J^j_)/, = (Jf3~J a +)fi- Complex conju- 
gation of (|A.6j) allows to find (J Q _ J^-)/;. 

When shells are the same, the result is as follows 

(J a+ J a+ ) fi = -Q] 1/2 (N + A)Q- 1/2 {N) J -^ exp(#a + iip' a ) tan(p' a /2) cot(p Q /2) 

N/2 

X ^(-) n (exp(#^ - iip a ) tan(/4/2) tan(p a /2)) n (iV + 2 - 2n)Qf l (N + 2-2n) (A.8) 

n=0 

= -QJ 1/2 (N + A)Q; 1/2 (N)j a (2j a + l)cxp(2^)tan 2 (p^/2) 

N/2 

x 2(-) n (exp(i#, - iip a )t<m(p'j2)t<m(p a /2)) n (n+ l)Q fi (N - 2n), (A.9) 

rc=0 

(J*+J«-) f i = ~Q] 1/2 (n)q; 1/2 (n)(^q^(n) 

N/2 

+ (j a + l/2)^(-)" +1 (cxp^/ - iip a )tan(p' a /2)tan(p a /2)) n (2nj a - l)Q fi {N - 2n)). (A.10) 



r 



The representations (|A.8|) and (|A.9[) are equivalent. a fixed number of the particles as well. However, the 

™ r „ . . number of shells is usually much larger than TV, and the 

1 he formalism allows also to provide an expression tor ,. u , , » r . , r , , 

^ 1 ratio between i I and N indicates the accuracy ol the so- 
the average as a sum over the quantum numbers a tor 
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lution of the variational problem. The higher the ratio, 
the higher the accuracy. For this reason, the summation 
over the number of particles is easier from the computa- 
tional point of view, and therefore preferable. Also, the 



recursive computation of the function Qfi(N) starting 
from small numbers of particles is more simple. 

The function Q fi (TV) can be written in the form [22[ 



QMN) = £ n { ^2-l a !mNj2)& ^^ + ^(P'J2) sinW2))^/ 2 

N! + N 2 + ...=N ot 

x (cos(p' Q /2) cos(p Q /2))°°/ 2 -^/ 2 . (A.ll) 



Here the summation is over all sets of occupation num- occupation numbers are even. In the special case of a 
bers N a , which gives the total number of particle N, the single shell, we reproduce the equation (|6.2j) . 
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